Análisis Causal de la Mortalidad de Renacuajos de Rana Carrizo

Métodos Analíticos

Autores/as

Blanca Estela García Manjarrez - 118886

Yuneri Pérez Arellano – 199813

Fecha de publicación

21 de mayo de 2025

Introducción

En este proyecto empleamos los datos experimentales de Vonesh & Bolker (2005) [1], quienes investigaron cómo la presencia de depredadores induce plasticidad en los tiempos de eclosión de la rana de junco africana Hyperolius spinigularis [2]. El estudio demostró que, aunque la eclosión anticipada reduce la mortalidad embrionaria hasta en un 40% bajo alta presión depredadora, las larvas resultantes —a pesar de su menor tamaño inicial— presentan una tasa de supervivencia incluso superior, gracias a respuestas compensatorias como aceleración del crecimiento y reducción de la fase acuática.

[1] Vonesh, J. R., & Bolker, B. M. (2005). Statistical tools for analyzing larval amphibian survival data. Ecology, 86(1), 172-182.

[2] Se refiere a la capacidad de los embriones de la rana para ajustar su desarrollo y eclosión en respuesta a cambios ambientales, como la presencia de depredadores o el secado de su hábitat

Los datos constan de 48 observaciones (“tanques de renacuajos”) con cinco variables clave:

Variable Descripción
density Densidad inicial de renacuajos
pred Factor indicador de presencia de depredadores
size Tamaño de los renacuajos
surv Número de renacuajos que sobrevivieron
propsurv Proporción de supervivencia (surv/density)

El conjunto de datos cuenta con información de 48 tanques clasificados en pequeños, medianos y grandes, dependiendo de la densidad de renacuajos en cada uno. Aunado a lo anterior, se dispone de datos sobre la supervivencia y la tasa de supervivencia en cada tanque.

Para nuestro análisis nos centraremos en surv como variable de respuesta (binomial) frente a density como total de ensayos, y consideraremos la estructura de agrupación en tanques como fuente de variabilidad no explicada por las manipulaciones experimentales.

Ignorar la heterogeneidad entre tanques (no pooling) o estimar un intercepto independiente para cada uno (efectos fijos) puede llevar, respectivamente, a subestimar o sobreestimar la variabilidad real. Por ello empleamos también un modelo bayesiano jerárquico (partial pooling) con interceptos variables, de este modo, cada tanque tiene su propio parámetro de línea base, y al mismo tiempo estimamos la dispersión entre tanques mediante un prior adaptativo, que aprende de los datos.

Nuestros objetivos son:

  1. Reproducir y extender los ejemplos de Statistical Rethinking aplicados a los datos de Reed Frogs.
  2. Modelar la mortalidad larval y las respuestas compensatorias posteriores usando modelos de efectos fijos (no pooling) y modelos multinivel (partial pooling).
  3. Evaluar la calidad y complejidad de cada modelo mediante diagnósticos MCMC y criterios de información (WAIC/LOO).
  4. Ilustrar el trade-off underfitting/overfitting mediante simulaciones con distintos tamaños de muestra.
  5. Desplegar un análisis causal formal con un DAG que recoja nuestros supuestos de identificación.

Con este enfoque buscamos profundizar en los costes y beneficios de la eclosión temprana inducida por depredadores, y demostrar cómo la regularización adaptativa de los modelos jerárquicos permite inferir efectos individuales de forma más robusta en presencia de datos jerarquizados y dispersos.

Datos

Estos datos muestran una gran variabilidad. Una parte obedece al tratamiento experimental, pero otra proviene de factores que no hemos medido. Imaginemos cada fila como un “tanque”, un microambiente donde viven renacuajos con características propias que no capturamos, y esas peculiaridades generan diferencias en la tasa de supervivencia, aun cuando las variables predictoras sean idénticas. En cada tanque realizamos múltiples observaciones (los renacuajos), por lo que tenemos medidas repetidas y heterogeneidad entre grupos.

Si obviamos la estructura de tanques y aplicamos un mismo intercepto para todos, podríamos pasar por alto variaciones esenciales en la supervivencia y enmascarar relaciones con otras variables. En cambio, asignar un intercepto distinto a cada tanque mediante variables dummy equivaldría a una “amnesia” constante: cada tanque es único, pero sus datos también nos ayudan a estimar el comportamiento de los demás. Olvidar por completo esa conexión, de un tanque a otro, no tiene sentido.

DAG

El modelo jerárquico propuesto busca reflejar las interacciones complejas entre la densidad, la supervivencia y otros factores fundamentales, con el fin de profundizar en cómo las ranas ajustan sus estrategias adaptativas frente al riesgo de depredación.

A continuación, se presenta un Diagrama Acíclico Dirigido (DAG) que ilustra las relaciones causales entre las variables de interés. Este DAG se basa en la premisa de que la densidad de renacuajos, el tamaño y la presencia de depredadores influyen en la supervivencia, y que la jerarquía de los tanques también afecta a estas relaciones.

DAG causal: densidad, tamaño, depredadores y jerarquía de tanque

Con

  • \(\textrm{T}=\textrm{Tanque}\)

  • \(\textrm{D}=\textrm{Densidad inicial}\)

  • \(\textrm{G}=\textrm{Tamaño}\)

  • \(\textrm{P}=\textrm{Depredadores}\)

  • \(\textrm{S}=\textrm{Supervivencia}\)

Aunque podríamos imaginar otras dependencias entre las variables, hay que tener presente que estos datos provienen de un experimento controlado. En un escenario natural, sería razonable investigar vínculos, como el efecto del tamaño de los renacuajos en la densidad poblacional, la influencia de los depredadores sobre esa densidad, o el papel de variables no registradas —por ejemplo, la disponibilidad de alimento u otros recursos— tanto en el tamaño como en la densidad, e incluso factores genéticos que modulen el desarrollo de los renacuajos. Sería muy valioso repetir estas estimaciones en cuerpos de agua naturales, en lugar de en tanques de laboratorio. Por ahora, este experimento nos permite centrarnos en la tasa de supervivencia bajo condiciones estrictamente controladas, estableciendo una base sólida para futuros estudios en la naturaleza.

Modelos

1. Modelo Total Pooling

En este primer modelo completamente agrupado asumimos que todos los tanques tienen la misma probabilidad de supervivencia.
No hay diferencias entre tanques, salvo la variación por la densidad inicial \(D_i\).

  • Datos
    • \(T\): número total de tanques (48).
    • Para cada tanque \(i\):
      • \(D_i\): densidad inicial (número de renacuajos)
      • \(S_i\): número de renacuajos supervivientes
  • Modelo \[ S_i \sim \textrm{Binomial}(D_i,p_i) \]

\[ \textrm{logit}(p_i) = \alpha_{T[i]} \]

\[ \bar{\alpha} = \textrm{Normal}(0, 1.5) \] \[ \alpha_j = \textrm{Normal}(\bar{\alpha}, 1.5) \]

Este modelo ignora la heterogeneidad entre tanques (total pooling) y servirá como línea base para comparar con el modelo jerárquico.

En el enfoque de complete pooling, asignamos una única \(\alpha\) a todos los tanques, de modo que el modelo asume idéntica probabilidad de supervivencia para los renacuajos en cada uno. Es como si tratáramos todos los tanques como réplicas exactas, sin captar ninguna heterogeneidad más allá de la variación provocada por la densidad inicial.

A continuación, mostramos los resultados obtenidos para la estimación de \(\alpha\) en nuestro modelo inicial.

parameter mean sd 2.5% 50% 97.5% n_eff Rhat
alpha 0.84 0.07 0.71 0.84 0.97 1,642.13 1.00
lp__ −685.69 0.73 −687.78 −685.41 −685.17 1,798.72 1.00

Y de la misma manera, podemos cotejar nuestras predicciones con los datos originales:

Predicciones del modelo 1 (total pooling)

Podemos ver que las predicciones (puntos rojos) siguen la media de supervivencia de los tanques (línea discontinua azul), lo que evidencia un subajuste del modelo. Los datos originales aparecen en gris.

2. Modelo No-Pooling

En este modelo no agrupado (no pooling) asignamos un intercepto \(\alpha_i\) distinto a cada tanque, pero no compartimos información entre ellos. Esto significa que cada tanque tiene su propio parámetro de línea base, y no hay aprendizaje entre ellos.

\[ S_i \sim \textrm{Binomial}(D_i,p_i) \]

\[ \textrm{logit}(p_i) = \alpha_{T[i]} \]

\[ \bar{\alpha} = \textrm{Normal}(0, 1.5)\] Lo cual se traduce en el siguiente código de Stan:

Podemos ver los resultados de nuestro segundo modelo para la estimación de \(\alpha\).

parameter mean sd 2.5% 50% 97.5% n_eff Rhat
alpha[1] 1.71 0.75 0.35 1.67 3.37 5,724.07 1.00
alpha[2] 2.40 0.90 0.82 2.35 4.31 5,328.56 1.00
alpha[3] 0.77 0.64 −0.42 0.77 2.03 5,472.50 1.00
alpha[4] 2.40 0.91 0.81 2.34 4.36 5,886.19 1.00
alpha[5] 1.72 0.75 0.33 1.67 3.27 5,334.79 1.00
alpha[6] 1.72 0.75 0.39 1.67 3.31 5,321.69 1.00
alpha[7] 2.42 0.94 0.81 2.37 4.46 5,569.04 1.00
alpha[8] 1.71 0.76 0.35 1.68 3.34 4,793.52 1.00
alpha[9] −0.38 0.62 −1.61 −0.36 0.78 6,062.49 1.00
alpha[10] 1.70 0.75 0.32 1.67 3.27 5,790.44 1.00
alpha[11] 0.78 0.65 −0.45 0.75 2.12 6,114.90 1.00
alpha[12] 0.37 0.61 −0.79 0.35 1.59 5,666.62 1.00
alpha[13] 0.76 0.62 −0.42 0.75 2.01 5,344.46 1.00
alpha[14] 0.00 0.59 −1.12 0.01 1.15 6,518.79 1.00
alpha[15] 1.70 0.75 0.35 1.66 3.28 5,225.61 1.00
alpha[16] 1.71 0.76 0.31 1.67 3.30 5,853.17 1.00
alpha[17] 2.54 0.67 1.37 2.49 4.01 5,191.84 1.00
alpha[18] 2.15 0.62 1.05 2.10 3.47 4,555.72 1.00
alpha[19] 1.80 0.54 0.86 1.78 2.95 4,589.60 1.00
alpha[20] 3.08 0.80 1.72 3.02 4.80 4,789.80 1.00
alpha[21] 2.15 0.59 1.10 2.10 3.40 4,753.37 1.00
alpha[22] 2.14 0.59 1.09 2.10 3.41 5,504.08 1.00
alpha[23] 2.15 0.59 1.10 2.12 3.38 4,821.24 1.00
alpha[24] 1.55 0.50 0.63 1.52 2.59 5,084.26 1.00
alpha[25] −1.10 0.46 −2.06 −1.09 −0.23 5,645.59 1.00
alpha[26] 0.08 0.40 −0.69 0.08 0.86 6,467.95 1.00
alpha[27] −1.54 0.51 −2.58 −1.51 −0.64 5,237.68 1.00
alpha[28] −0.55 0.41 −1.37 −0.54 0.22 5,866.64 1.00
alpha[29] 0.07 0.37 −0.68 0.07 0.80 5,883.37 1.00
alpha[30] 1.30 0.47 0.40 1.28 2.24 6,271.42 1.00
alpha[31] −0.72 0.43 −1.56 −0.72 0.10 5,240.80 1.00
alpha[32] −0.39 0.41 −1.21 −0.38 0.39 5,668.79 1.00
alpha[33] 2.84 0.67 1.67 2.80 4.24 5,775.01 1.00
alpha[34] 2.46 0.60 1.43 2.42 3.77 5,478.55 1.00
alpha[35] 2.48 0.59 1.43 2.45 3.76 4,548.57 1.00
alpha[36] 1.91 0.47 1.05 1.88 2.91 5,854.25 1.00
alpha[37] 1.90 0.47 1.06 1.88 2.90 5,773.73 1.00
alpha[38] 3.37 0.79 2.00 3.30 5.13 5,476.40 1.00
alpha[39] 2.48 0.59 1.44 2.44 3.70 4,876.90 1.00
alpha[40] 2.16 0.52 1.21 2.14 3.29 4,951.47 1.00
alpha[41] −1.91 0.49 −2.95 −1.88 −1.01 5,298.03 1.00
alpha[42] −0.64 0.35 −1.36 −0.64 0.03 5,494.54 1.00
alpha[43] −0.51 0.35 −1.20 −0.50 0.16 4,819.55 1.00
alpha[44] −0.39 0.33 −1.04 −0.38 0.25 5,822.05 1.00
alpha[45] 0.52 0.33 −0.12 0.51 1.20 6,133.78 1.00
alpha[46] −0.64 0.34 −1.33 −0.63 0.01 5,374.62 1.00
alpha[47] 1.91 0.49 1.04 1.89 2.96 5,379.99 1.00
alpha[48] −0.06 0.34 −0.71 −0.05 0.61 6,345.95 1.00
lp__ −524.49 5.15 −535.51 −524.15 −515.51 1,780.74 1.00

Predicciones del modelo 2 (no pooling)

Al evaluar las estimaciones del modelo sin pooling (un intercepto \(\alpha_i\) independiente por tanque) obtenemos:

  1. Gran variabilidad entre tanques
    • Los \(\alpha_i\) (en escala log-odds) oscilan aproximadamente entre \(-1.9\) y \(+3.4\).
    • Transformados a probabilidad, algunos tanques se estiman con supervivencias cercanas al 10–20% y otros al 90–95%.
  2. Incertidumbre muy desigual
    • Tanques con pocas observaciones (densidad pequeña) presentan desviaciones estándar de \(\alpha_i\) de 0.6–0.8 y rangos de credibilidad muy amplios (p. ej. \(\alpha_3: [–0.41, +2.06]\).
    • Tanques con más renacuajos reducen su incertidumbre a 0.3–0.5 en la desviación estándar de \(\alpha_i\).
  3. Sobreajuste
    • Las predicciones del modelo (puntos rojos) siguen casi exactamente los datos observados (puntos grises), incluso en valores extremos.
    • No existe “arrastre” hacia un promedio general: cada tanque se ajusta únicamente con su propia información.
  4. Problemas en tanques pequeños
    • Con muestras muy pequeñas, pocas muertes o supervivencias cambian drásticamente la estimación de \(\alpha_i\).
    • El ancho de los intervalos de credibilidad hace poco útiles esas predicciones para la toma de decisiones.

Por lo anterior, el modelo no agrupado captura fielmente cada dato empírico, pero padece de sobreajuste y de alta incertidumbre en tanques con pocas observaciones. Para obtener estimaciones más estables y evitar extremos sin fundamento, es recomendable utilizar un modelo jerárquico (partial pooling) que comparta información entre tanques.

3. Modelo Partial Pooling

En este tercer modelo parcialmente agrupado (partial pooling) asignamos un intercepto distinto a cada tanque, pero también estimamos la variabilidad entre ellos. Esto nos permite captar la heterogeneidad entre tanques y, al mismo tiempo, compartir información entre ellos.

\[ S_i \sim \textrm{Binomial}(D_i,p_i) \]

\[ \textrm{logit}(p_i) = \alpha_{T[i]} \]

\[ \bar{\alpha} = \textrm{Normal}(0, 1.5)\]

\[ \sigma_\alpha = \textrm{Exponential}(1) \]

\[ \alpha_j = \textrm{Normal}(\bar{\alpha}, \sigma_\alpha) \]

Lo cual se traduce en el siguiente código de Stan:

parameter mean sd 2.5% 50% 97.5% n_eff Rhat
alpha_tank[1] 2.14 0.87 0.59 2.08 4.01 4,295.20 1.00
alpha_tank[2] 3.08 1.12 1.17 2.97 5.52 4,108.25 1.00
alpha_tank[3] 1.00 0.67 −0.24 0.97 2.43 5,376.00 1.00
alpha_tank[4] 3.08 1.12 1.20 2.97 5.66 3,948.65 1.00
alpha_tank[5] 2.13 0.87 0.63 2.07 4.02 5,107.69 1.00
alpha_tank[6] 2.14 0.88 0.64 2.08 4.01 4,268.95 1.00
alpha_tank[7] 3.07 1.09 1.29 2.97 5.52 4,277.00 1.00
alpha_tank[8] 2.13 0.87 0.64 2.08 4.01 4,481.12 1.00
alpha_tank[9] −0.16 0.60 −1.34 −0.15 0.98 5,136.20 1.00
alpha_tank[10] 2.13 0.86 0.61 2.07 3.99 4,134.62 1.00
alpha_tank[11] 1.01 0.66 −0.24 0.98 2.33 5,138.46 1.00
alpha_tank[12] 0.57 0.63 −0.61 0.56 1.87 5,793.49 1.00
alpha_tank[13] 1.00 0.66 −0.21 0.97 2.40 5,000.17 1.00
alpha_tank[14] 0.19 0.60 −0.97 0.19 1.39 4,957.28 1.00
alpha_tank[15] 2.15 0.89 0.57 2.09 4.06 4,268.56 1.00
alpha_tank[16] 2.14 0.87 0.63 2.07 4.06 4,682.98 1.00
alpha_tank[17] 2.91 0.77 1.58 2.84 4.61 4,397.65 1.00
alpha_tank[18] 2.39 0.68 1.20 2.34 3.88 5,630.03 1.00
alpha_tank[19] 2.01 0.58 0.99 1.96 3.25 4,262.17 1.00
alpha_tank[20] 3.68 1.04 1.98 3.56 6.19 3,480.68 1.00
alpha_tank[21] 2.40 0.65 1.22 2.36 3.80 3,978.77 1.00
alpha_tank[22] 2.39 0.64 1.25 2.36 3.75 4,613.55 1.00
alpha_tank[23] 2.40 0.66 1.23 2.36 3.85 4,557.86 1.00
alpha_tank[24] 1.70 0.53 0.76 1.67 2.81 4,068.18 1.00
alpha_tank[25] −0.99 0.44 −1.91 −0.98 −0.18 4,922.89 1.00
alpha_tank[26] 0.17 0.40 −0.62 0.17 0.97 3,935.60 1.00
alpha_tank[27] −1.43 0.50 −2.47 −1.41 −0.50 4,518.75 1.00
alpha_tank[28] −0.47 0.41 −1.31 −0.48 0.35 4,989.20 1.00
alpha_tank[29] 0.16 0.39 −0.59 0.16 0.93 5,348.16 1.00
alpha_tank[30] 1.44 0.50 0.53 1.43 2.47 5,257.82 1.00
alpha_tank[31] −0.63 0.41 −1.45 −0.62 0.15 4,954.26 1.00
alpha_tank[32] −0.30 0.40 −1.09 −0.29 0.45 5,028.63 1.00
alpha_tank[33] 3.18 0.73 1.91 3.11 4.79 4,527.59 1.00
alpha_tank[34] 2.71 0.64 1.60 2.66 4.14 3,872.89 1.00
alpha_tank[35] 2.72 0.64 1.60 2.68 4.14 3,791.94 1.00
alpha_tank[36] 2.08 0.52 1.13 2.05 3.16 4,543.74 1.00
alpha_tank[37] 2.06 0.50 1.16 2.02 3.14 4,933.04 1.00
alpha_tank[38] 3.90 0.96 2.32 3.79 6.07 3,566.87 1.00
alpha_tank[39] 2.70 0.63 1.63 2.64 4.11 3,776.54 1.00
alpha_tank[40] 2.35 0.58 1.30 2.31 3.60 4,673.15 1.00
alpha_tank[41] −1.80 0.47 −2.81 −1.78 −0.95 5,278.21 1.00
alpha_tank[42] −0.58 0.35 −1.28 −0.57 0.10 4,986.06 1.00
alpha_tank[43] −0.45 0.35 −1.17 −0.45 0.20 5,666.27 1.00
alpha_tank[44] −0.33 0.34 −1.02 −0.33 0.33 5,049.41 1.00
alpha_tank[45] 0.58 0.36 −0.11 0.57 1.29 5,148.54 1.00
alpha_tank[46] −0.57 0.34 −1.24 −0.55 0.09 5,300.03 1.00
alpha_tank[47] 2.08 0.52 1.17 2.04 3.18 4,949.98 1.00
alpha_tank[48] 0.00 0.33 −0.65 0.00 0.67 5,287.28 1.00
lp__ −532.64 5.47 −544.33 −532.31 −522.85 1,300.14 1.00

Y para los hiperparámetros

parameter mean sd 2.5% 50% 97.5% n_eff Rhat
mu_alpha 1.35 0.26 0.86 1.34 1.87 3,344.76 1.00
sigma_alpha 1.62 0.21 1.25 1.60 2.09 2,729.90 1.00

Predicciones del modelo 3 (partial pooling)
  1. Hiperparámetros poblacionales
    • Media global \(\mu_\alpha\approx1.35\) (sd = 0.25 - supervivencia promedio) \(\text{logit}^{-1}(1.35)\approx0.79\).
    • Variabilidad entre tanques \(\sigma_\alpha\approx1.63\) (sd = 0.22), muestra heterogeneidad moderada-alta en tasas de supervivencia.
  2. Regularización adaptativa (“shrinkage”)
    • Las estimaciones por tanque (\(\alpha_j\)) se encogen hacia la media global:
      • Más encogimiento en tanques pequeños (poca información).
      • Menos encogimiento en tanques grandes (más datos).
  3. Equilibrio under-/over-fitting
    • El modelo aprende el grado óptimo de pooling desde los datos.
    • Corrige el subajuste del complete-pooling y el sobreajuste del no-pooling.
  4. Precisión y convergencia
    • Intervalos más estrechos que en el modelo no-pooling, pero más amplios que en el total-pooling.
    • Todos los \(\hat R\approx1\) y n_eff elevados garantizan buena convergencia de la cadena.
  5. Implicaciones ecológicas
    • Capta diferencias reales entre tanques (densidad, depredadores), sin ignorar heterogeneidad.
    • Cada tanque “aprende” de los demás, mitigando el sesgo de muestras pequeñas.

El modelo parcialmente agrupado ofrece el mejor compromiso: reduce el ruido de tanques con pocos datos y a la vez preserva las verdaderas diferencias entre ellos, obteniendo estimaciones más robustas y ecológicamente interpretables.

Comparativos

Diagnósticos MCMC

WAIC / LOO

Predicciones por tanque

Simulación y Validación (Trade-off Under/Overfitting)

1. Generar datos simulados

2. Ajustar modelo jerárquico (partial pooling)

3. Gráfico de errores (Figura 12.3)

Conclusiones